home *** CD-ROM | disk | FTP | other *** search
/ SPACE 1 / SPACE - Library 1 - Volume 1.iso / program / 441 / fplib20 / itrig.c < prev    next >
C/C++ Source or Header  |  1990-11-23  |  2KB  |  98 lines
  1. /* ASIN, ACOS, ATAN, and ATAN2 functions for Sozobon C.            */
  2. /* Copyright © David Brooks, 1989 All Rights Reserved            */
  3.  
  4. #include <fplib.h>
  5. #include <errno.h>
  6.  
  7. /* ARCTAN:                                */
  8. /*                                    */
  9. /* The method is from Hart et al: "Computer Approximations".  It uses a    */
  10. /* divide-and-conquer algorithm.  The argument range is divided into    */
  11. /* five using the partition points tan({1,3,5,7}*pi/16), a polynomial    */
  12. /* valid over +/- pi/16 is calculated, and other magic is used to    */
  13. /* reposition the result.  Precision is >8 places.            */
  14.  
  15. float atan(a) fstruct a;
  16. {    fstruct        absval;
  17.     register float    tx0, tsq, temp;
  18.     register int    part;
  19.     register char    sign;
  20.  
  21.     static float    adj[] = {0.0, 0.414213562, 1.0, 2.414213562};
  22.     static float    atof[] = {0.0, 0.39269908, 0.78539816,
  23.                     1.178097245, 1.57079633};
  24.  
  25.     sign = a.sc[3] & 0x80;
  26.     a.sc[3] &= 0x7f;            /* get fabs(a) */
  27.     tx0 = a.f;
  28.  
  29. /* Figure out the partition */
  30.  
  31.     part = tx0<0.66817864?(tx0<0.19891237?0:1):(tx0<1.49660576?2:
  32.             tx0<5.0273395?3:4);
  33.     if (part == 4)
  34.         tx0 = -1.0 / tx0;
  35.     else if (part != 0)
  36.     {    temp = adj[part];
  37.         tx0 = (tx0 - temp) / (tx0 * temp + 1.0);
  38.     }
  39.  
  40. /* Here's the calculation */
  41.  
  42.     tsq = tx0 * tx0 + 1.67784279;
  43.     a.f = (0.93833093 / tsq + 0.44075154) * tx0 + atof[part];
  44.  
  45.     if (a.sc[3] != 0)            /* Negate if negative */
  46.         a.sc[3] ^= sign;
  47.     return a.f;
  48. }
  49.  
  50. /* ARCSIN and ARCCOS use the standard identities.  There's some less    */
  51. /* optimal code here becausethe released Sozobon C can't properly    */
  52. /* compare two negative floating numbers.                */
  53.  
  54. float asin(a) float a;
  55. {
  56.     if (a > 1.0 || (a+1.0) < 0.0)
  57.     {    errno = EDOM;
  58.         return 0.0;
  59.     }
  60.  
  61.     return atan(a / sqrt(1.0 - a*a));
  62. }
  63.  
  64. float acos(a) float a;
  65. {    float    temp;
  66.  
  67.     if (a > 1.0 || (a+1.0) < 0.0)
  68.     {    errno = EDOM;
  69.         return 0.0;
  70.     }
  71.  
  72.     temp = atan(sqrt(1.0 - a*a) / a);
  73.     if (a >= 0.0)
  74.         return temp;
  75.     else
  76.         return M_PI + temp;
  77. }
  78.  
  79. /* ATAN2:                                */
  80. /*                                    */
  81. /* Computes atan(quotient), and returns that for positive cos, else     */
  82. /* extends the range depending on the sin.                */
  83.  
  84. float atan2(s, c) fstruct s, c;
  85. {    register float    r;
  86.  
  87.     if (c.sc[3] == 0)            /* Infinite argument */
  88.         return (s.sc[3]<0)?-M_PI_2:M_PI_2;
  89.  
  90.     r = atan(s.f / c.f);
  91.  
  92.     if (c.sc[3] >= 0)
  93.         return r;            /* Range -pi/2..pi/2 */
  94.     if (s.sc[3] >= 0)
  95.         return M_PI + r;
  96.     return -M_PI + r;
  97. }
  98.